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Abstract. In this work, we present a conditionally stable finite-difference 
scheme that consistently approximates the solution of a general class of (3 + 1)- 
dimensional nonlinear equations that generalizes in various ways the quantita- 
tive model governing discrete arrays consisting of coupled harmonic oscillators. 
Associated with this method, there exists a discrete scheme of energy that con- 
sistently approximates its continuous counterpart. The method has the prop- 
erties that the associated rate of change of the discrete energy consistently 
approximates its continuous counterpart, and it approximates both a fully 
continuous medium and a spatially discretized system. Conditional stability 
of the numerical technique is established, and applications are provided to the 
existence of the process of nonlinear supratransmission in generalized Klein- 
Gordon systems and the propagation of binary signals in semi-unbounded, 
three-dimensional arrays of harmonic oscillators coupled through springs and 
perturbed harmonically at the boundaries, where the basic model is a modified 
sine-Gordon equation; our results show that a perfect transmission is achieved 
via the modulation of the driving amplitude at the boundary. Additionally, 
we present an example of a nonlinear system with a forbidden band-gap which 
does not present supratransmission, thus establishing that the existence of a 
forbidden band-gap in the linear dispersion relation of a nonlinear system is 
not a sufficient condition for the system to present supratransmission. 



1. Introduction 

Almost five years after the appearance of the pioneering letter by Geniet and 
Leon PP, the phenomenon of nonlinear supratransmission has been studied widely 
in many one-dimensional, physical systems. The phenomenon consists in a sudden 
increase in the energy injected into a nonlinear system by a harmonic perturbation 
irradiating at a frequency in the forbidden band-gap, and the research in the field 
has concentrated mainly on discrete media such as mechanical chains of oscillators 
described by coupled sine-Gordon and Klein-Gordon equations pQ, coupled double 
sine-Gordon equations [2J, Fermi-Pasta-Ulam nonlinear chains [3J, and Bragg media 
in the nonlinear Kerr regime [4]. Nonetheless, some research has been done in 
the continuous case scenario, where the sine-Gordon equation has been a common 
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denominator [5l|6]. Meanwhile, from a pragmatic perspective, the importance of 
the process of nonlinear supratransmision has been evidenced through the many 
applications proposed to the design of digital amplifiers of ultra weak signals [7], 
light detectors sensitive to very weak excitations jS], optical waveguide arrays [9], 
and light filters [10] . 

Of course, the problem in the numerical study of the process of nonlinear supra- 
transmission lies in the development of a reliable computational technique to ap- 
proximate consistently the solutions to the mixed-value problem, the local energy 
density of the system, and its total energy, in view of the fact that supratransmis- 
sion is better characterized in the energy domain. Moreover, from a historically 
point of view the use of symplectic methods for Hamiltonian systems has proved 
to yield more than satisfactory results |llj : unfortunately, the medium we analyze 
in the present work contemplates the inclusion of internal and external damping 
terms which make it nonconservative in general. Nonetheless, the main part of our 
study will be devoted to develop a finite-difference scheme for the problem under 
analysis, together with a discrete scheme for the local energy density and the total 
energy of the system with consistency properties not only in the energy domain, 
but also in the domain of the rate of change of energy of the medium. 

In general, the study of (3 + l)-dimensional systems governed by sine-Gordon 
equations is an important problem in the physical sciences. For instance, ring- 
shaped solitary wave solutions of these type of systems have been numerically in- 
vestigated to show ultimately that such solutions have quasi-soliton properties [12] . 
The existence of multi-soliton and vortex-soliton solutions has been established for 
this model, too [13]. N- layer sine- Gordon- type models have been studied in or- 
der to generalize the results obtained for the two-layer sine-Gordon model [Tl], a 
model that has been used to describe the dynamics of high transition temperature 
superconductors [15) . Finally, the (3 + l)-dimensional sine-Gordon equation has 
been used to explore the possibility of stable superluminal propagation of short 
electromagnetic excitations [16]. 

In Section [5] of this work, we present the (3 + l)-dimensional problem under 
study in its most general form. The model includes the presence of internal and 
external damping, relativistic mass, and generalized Josephson currents. Here, 
we present the Lagrangian of the undamped case as well as an energy analysis 
of the system under study, and a statement of a similar problem in spherically 
symmetric media. Section [3] introduces the finite-difference schemes employed to 
approximate solutions of the mixed- value problem of interest, and the schemes used 
to approximate the local energy density and the total energy of the system. We 
establish that the discrete rate of change of energy is a consistent estimate of its 
continuous counterpart, and a stability condition is proved. In Section |4] we show 
numerically that the process of nonlinear supratransmission is present in the semi- 
discrete system under scrutiny by means of an application of the method presented 
in this work. The relevance of our results will be shown when we demonstrate next 
that the system under study does not support supratransmission when the medium 
is radially perturbed at the origin, whence it will follow the existence of a forbidden 
band-gap for the frequency in the linear dispersion relation of a nonlinear system 
does not necessarily guarantee the presence of supratransmission in the medium. A 
second application to the generation and propagation of localized nonlinear modes 
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in the system of interest is presented next, and we close our work with a section of 
concluding remarks. 



In the present section we introduce the two mathematical models under study in 
this work. Throughout, the nonnegative constants /3 and 7 represent, respectively, 
the coefficients of internal and external damping of the medium, and the pure-real 
or pure-imaginary constant m denotes a relativistic mass; this last parameter has 
been included to suggest further applications of our results to the field of particle 
physics |17) . Moreover, the nonnegative value J will be called generalized Josephson 
current, and its inception has been realized with applications to superconductivity 
in mind [18] . 

2.1. Cartesian problem. Let us represent the closure of the first octant of the 
Euclidean space M 3 by D, let V be any continuously diffcrcntiable real function 
defined in all of R, and assume that u is a function of (x, t), where x € D and 
t € R + . Under these circumstances, the medium studied in the present paper is 
described by the generalized partial differential equation 



in which V 2 represents the Laplacian operator. 

It is important to point out that ([lj generalizes nonlinear partial differential 
equations such as the Klein-Gordon equation, the sine-Gordon equation, and the 
Landau-Ginzburg equation, amongst others. If m is a pure-real number then a 
modified sine-Gordon model is obtained, for instance, when a potential of the form 
V(u) = 1 — cos it is considered, and a modified nonlinear Klein-Gordon equation 
results when V(u) = ^u 2 — -gu 4 . Meanwhile, for every positive real number A, a 
modified Landau-Ginzburg equation is obtained if V(u) — Ait 4 when m is a pure- 
imaginary number. 

This investigation considers particularly the study of the sine-Gordon and Klein- 
Gordon equations, two models that have been thoroughly studied in the literature 
PH HOI HD H3 123 HI- The inclusion of the parameters /3 and 7 in our model 
correspond to the need of considering generalizations of physically realistic models 
in which internal and external damping are present, such as problems arising in the 
study of long Josephson junctions between superconductors when dissipative effects 
are taken into account [18] or in the investigation of fluxons in Josephson trans- 
mission lines [25] ■ Mathematically, the study of sine-Gordon and Klein-Gordon 
systems where linear damping is present has lead to the discovery of weak solu- 
tions of these equations [26], the proof of the existence of the maximal attractor 
in dissipative systems of Klein-Gordon-Schrodinger equations [27] , the discovery of 
the mechanism of the ratchet-like dynamics of solitons in dissipative Klein-Gordon 
media driven by a bi-harmonic force |28j . the proof of the existence of multista- 
bilities and soliton trapping in the damped Klein-Gordon equation with external 
periodic excitation via the asymptotic perturbation method [29], amongst many 
other analytical results [30 ] I3T ] 132 ] 133 ] 131 ] . 

In the case of conservative sine- Gordon and Klein- Gordon media described by ([1]) 
with a generalized Josephson current equal to zero, the linear dispersion relation is 
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obtained when considering solutions in the linearized systems in the form of linear 
modes 35 . In such cases, the dispersion relation adopts the form 



(2) o.'-ii.^.//) m 2 + 1 + 4 ( sin 2 f +sin 2 | +sin 2 V 




which possesses a forbidden band-gap given by f2 < \/m 2 + 1. From a practical 
point of view, the parameter ft will represent the frequency of the driving boundary 
in a harmonically perturbed system described by (H}. 

Once the pragmatic importance of our model has been understood, we proceed 
to simplify it by letting 

(3) G(u) = ^rn 2 u 2 + V{u) - Ju. 

In these terms, the Lagrangian associated with the conservative portion of ([I]) and 
the corresponding Hamiltonian are 

(4) 

-G(u) and n = -i^—j + -\\\7u\\ 2 + G(u), 

respectively, where || • || represents the Euclidean norm in R 3 . Moreover, it is easy 
to derive the following expression for the total energy of the system at any time t, 
in which the integrand is the local energy density: 

(5) E(t) = JJJ Hdx. 

Here, it is important to observe that the total energy of our problem is positive 
whenever G is a nonnegative real function, for instance, when the system has no 
generalized Josephson current and V(u) — u p with p an even positive integer. 

For computational reasons, we will restrict our study to bounded domains D of 
the form [0,L] x [0,L] x [0,L], where L is a positive constant. Moreover, we will 
assume that Neumann boundary data will be imposed on the sides of D opposite 
to the origin. More precisely, we will assume that Vu ■ n = on the sides of D 
opposite to the origin. Furthermore, it will be important to consider Dirichlet data 
of the form u(x, t) = 4>{t), where x are on the sides of D adjacent to the origin. 

Proposition 1. The instantaneous rate of change with respect to time of the total 
energy associated with the partial differential eguation {1} in the region D = [0, L] 3 
of R 3 with boundary data Vu ■ n = on the sides of D opposite to the origin and 
Dirichlet condition on the sides dD + adjacent to the origin, is given by 



(6) 









IL 




dxj 









NUMERICAL METHOD WITH ENERGY CONSISTENCY FOR SINE-GORDON 



5 



Proof. Taking derivative with respect to time on both sides of ([5]), using Green's 
first identity, and substituting equation ([T]) next, we obtain that 



E'{t) 



du 
~dt 



d 2 u 



+ G'{u)}dx + 



D 



dt 



Vu\\ 2 dx 



du ( d 2 u 



dt dt 2 



P 



D 



dt 



- V 2 u + G'(u) }dx 



du 

It. 



Vu • n da 



dt 



dD 



du 

— Vm • n da. 

3D Ot 



On the other hand, from Green's first identity we see that 



D 



du 
~di 



du 
l)t 



dx = 



du d 



(Vu • fi) da — 



du 

~dt 



dx. 



idD dt dt 

The surface integrals in these last two equations are equal to zero on the three sides 
of D opposite to the origin, whence the result follows. □ 

It is worth noticing that, in view of the hypotheses of Proposition [1] more con- 
crete expressions for some terms in ([6]) are readily at hand. Particularly, it is 
convenient to observe that 



0D + 



du 
It 



Vm ■ iida = — 



(7) 



L du 9ti 
-g^{0,y,z) — (0,y,z)dydz 

157 (a, 0, z)^(x, 0, z) dx dz 
dt dy 

~fo( x ' y> °)^( ;E ' y ' °) dxd V- 



Similarly, 
(8) 



du d „. 

, w+ md-t iVu - n)da 



r i 9u, n . 



dt dx 
d 2 u 



dtdy 



(0,y,z)dy dz 
(x, 0, z) dx dz 
(x,y,0) dx dy. 



dt dz 



It is also important to notice that if j3 = 7 = and if either 



dt 



or Vm ■ h = 



on dD + , then the energy of the system is conserved throughout time. 

2.2. Spherical problem. As usual, let D be the closure of the first octant in M 3 , 
and assume that u is a radially symmetric solution of (TTJ. Let r = ||x|| represent 
the Euclidean norm of the vector x 6 R 3 , let v(r,t) = ru(r,t), and assume that 
u is a solution of problem (TTJ) for Dirichlet boundary data in the origin given by 
u(0,t) — 4>(t). Then v satisfies the relation v(0,t) = 0, together with the partial 
differential equation 

, , d 2 v d 2 v 9 , d 3 v dv 

Computationally and for the remainder of the present section, the region D will 
represent the closure of the portion of the solid sphere with center in the origin and 
radius equal to L that lies in the first octant, and Neumann boundary data will 
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be imposed on the boundary of the region. Moreover, since the Dirichlet boundary 
condition of u at the origin translates into a void condition for v (which in turn 
translates into a trivial solution for problem ^ when vanishing initial conditions 
are chosen and the Josephson current is equal to zero), we set f(e,i) = t(f>{t) for 
some e > sufficiently close to zero. 

Under the presence of spherical symmetry and assuming that u{L 1 1) = for 
every t € M + , the energy expression ([5]) of the undamped system may be computed 
in terms of v via the following expression: 




This follows immediately after noticing that (|^) = (f*§^) + -§^{ru 2 ). Here, G 
is the simplified form of the potential function provided in the previous subsection. 
Moreover, in view of Proposition [TJ a rate of change of the total energy of system 
(O is readily established as our next analytical result. 

Proposition 2. The instantaneous rate of change of the total energy of a system 
satisfying Q on the set D, with Dirichlet boundary condition at the origin, and 
Neumann data Vw • h = and Dirichlet data u(L,t) — in the intersection of D 
with the set o/x£l 3 such that ||x|| = L, is provided by the formula 




Proof. It follows directly from Proposition [T] and the substitution v(r,t) — ru(r,t). 

□ 



In this case we must observe that the Neumann boundary data will take the 
form ^ = on the curved side of the wedge D. In terms of the variable v, this 
condition translates into the equation 

(12) ^ + -=0. 

or r 



2.3. Discrete problem. In this section, we introduce a model that describes the 
dynamics of a discrete system of pendula attached springs. Let w m! „ iP be a real 
function on the real variable t, for every m,n,p € Z + U {0} and every t > 0. We 
will consider now the infinite system of coupled ordinary differential equations with 
constant coupling coefficient c > 0, in which m, n,p £ Z + : 

(13) Um : n,p c \7 u m ^ n ^p ~t- m u m ^ n p -(- V (u m , nj p) J — (3^\7 ii ra n p 7^,fi,p 

Here, the discrete Laplacian operator V 2 , the Hamiltonian H of the lattice site at 
position (to, n,p) for the conservative theory, and the total energy E of the system 
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are, respectively, 



' N N 

{ u l,n,p — u Q,n,p) 2 + J]] ( u m,l,p ~ U-mfitf) 
n^p—1 m,p= 1 



-^m,n,p — 7^ {^m,n.p ^ (^m+l,n,p ^m,n,p) ~t~ C (^m,n+l,p ^m,n,p) 

(14) -j-C (tXm,n,p+l ^m,n,p) rTl ^ m5n; p} "i" ^(^m.n,p) <-^m,n,p: 

N 

E ^ ^ Hrn.n,p ~\~ 

m,n,p— 1 
AT 

m,n— 1 

2 

It is important to remark here that the inclusion of the terms multiplied by 
in the discrete energy corresponds to the need to include the potential from the 
coupling of the nodes adjacent to the boundary. 

For the sake of convenience, we introduce the notation 

^x^m,n,p — ^m+l,n,p ^m.n,p: 
(15) ^y^m,n.p — ^m,n+l,p ^m,n,p> 



&z^m.n,p ^m,n,p+l ^m,n,p- 

Moreover, for computational reasons we will assume that m, n and p take on values 
in the set {0, 1, . . . , N + 1} for a relatively large positive integer N, and assume 
that discrete Neumann boundary data are imposed on the boundaries n = N + 1, 
m = N + 1 and p = N + 1, that is, we assume that 

(16) 5 x UN,m,p — &yUm,N,p = ^zU m ,n,N = 0, 

for every m,n,p G {1, . . . , N}. Meanwhile, Dirichlet data will be required on the 
remaining boundaries. 

System (Tl3|) describes the evolution of a semi-unbounded, three-dimensional ar- 
ray of harmonic oscillators coupled through identical springs with a coupling co- 
efficient equal to c. The pendula are located at the discrete sites (r7i,n,p), where 
m,n,p € Z + , and the springs are parallel to a coordinate axis. Evidently, site 
(m,n,p) is coupled with the six sites (m ± l,n,p), (m,n± l,p) and (m,n,p± 1), 
and c represents the common coupling coefficient. Moreover, in practice we will 
subject the oscillators on the boundaries to harmonic driving in the form of the 
Dirichlet conditions 

(17) W m ,n,0 = Um,0,p = «0,n,p = Asm(fM), 

where Q is a frequency in the forbidden band-gap of the continuous-limit medium. 
A schematic representation of such a system is depicted in Fig. [T] In this context, it 

2 

is important to notice that the term with coefficient ^- in the expression E for the 
total energy of the system corresponds to the potential energy due to the coupling 
to the driving boundary. 

Our next result is a one-dimensional version of Green's first identity. For a proof, 
we refer to 136 . 
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Figure 1 . Schematic representation of a medium governed by sys- 
tem (fT5)) . The nodes located in sites (m, n,p) for m, n,p E Z + U{0} 
may physically represent harmonic oscillators, and the attaching 
segments of line play the roles of springs. 



Lemma 3 (Green's discrete first identity). For every sequence (a n ) n =Q for which 
a-N+i — cln, 



N 



N 



/Xa n +i ~ 2a„ + a„_i)a„ = a (a - ai) - /,(«r» ~ a n-i) 



□ 



Proposition 4. Consider a system satisfying (|f 3[) form, n,p = 1,2, . . ., iV+1, sub- 
ject to discrete Neumann conditions of the form b~ x iiN,m,p = $yU m ,N,p — ^z u m,n,N 
on the boundaries n = JV+1, m = iV + 1, andp = N + 1, and subject to Dirichlet 
data on the remaining boundaries. Then, the instantaneous rate of change of the 
energy of the node in site (m,n,p) with respect to time is given by 



dE 



N 



N 



lm,n,p-l 



A 



+ X] [(d x u Q ,ij)u ,ij + (S y -Uifl,j)ui,o,j + (<J*«i,j,o)«t,i,o] > - 7 ^ ( 

i,j=l J m,n,p— 1 
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" (Wm+l,n,p ^rn,n,p) — (^m,n,p -^m+l,n,p) C (lij7i-j-l,n,p 2lZ mjn .pH-?i m _ l,n.p)'U»m,n,p- 



Proof. Define Km,n,p — C ^rn,n,p{ < f U'rn,n,p ^m— l,n,p)) -^-m,n.p — ^ Urn 7 n.p(Um,n,p 

V»-i,p) and K m,n, P = ~c 2 u m , n ,p{u m ,n,p-u m ,n, P -i), for every m,n,pe{l,2,..., TV}. 
It is necessary to observe first of all that 

2 dt' 

Similar relations may be obtained for the derivatives of the other terms in the 
Hamiltonian which are multiplied by c 2 . Moreover, taking derivative of the Hamil- 
tonian with respect to t, one obtains that 

j£ = {K m ,n,p ~ ^m+l,n,p) + (K m ,n,p ~ ^m,n+l,p) + (^m,n,p ~ ^m,n,p+l) 

~t~ {^m,n,p C V7 Um,n,p ~t~ TTX Um.n,p ~\~ V (Um,n.p) J } Um,n.p 
= {K m ^ lt p — -ftT m +l : n,p) + {K m n p — -^ m ^ n +l !?) ) + {K m n p ~~ K m n p+1 ) 

~\~fi\ilrn J c\,n,p '^'U"m,n,p ~t~ ^m— l,n,p)^m,n,p ~t~ j^{^rn,n+l.p ^^rn.n,p ^m.n— l,p)^m,n,p 

and sum over indexes of m, n and p in the set {1, 2, . . . , iV}. We identify the sums 
of the first three expressions in parenthesis as telescoping series and proceed to 
simplify; at the same time, three applications of the discrete version of Green's first 
identity provide alternative expressions for the terms multiplied by (3. On the other 
hand, by differentiating the energy expression with respect to time, substituting the 
derivative of the Hamiltonians and simplifying, we reach the desired formula. □ 



3. Numerical analysis 



3.1. Cartesian problem . In order to approximate solutions of the partial differ- 
ential equation ^ on the cube [0, L] x [0, L]x [0, L] over an interval of time of length 
T, we choose a regular partition = to < ti < • ■ • < tu = T of [0, T] with time step 
equal to At, as well as three regular partitions of [0, L] consisting of N x + 1, N y + 1 
and N z + 1 subintervals, each with step equal to Ax, Ay and Az, respectively. For 
all permissible indexes k, m, n and p, we represent the approximate solution to 
our problem at time kAt and at the location (mAi, nAy,pAz) by u: 



m 7 n,p' 



The dis- 
cretization of the problem under study is provided by the finite-difference schemes 
(19) 



6?ut 



(At) 2 



- 1 



(Ax) 2 (Ay) 2 (Az) 2 



m 

~2~ 



,k+l 



k-l 



m,n,p ' 2At m ' n 'P 

v(u k+1 ) - v(u k -l ) 



t^m,n.p 



.k-l 



for every k = 1, . . 
to the conditions 



, M — 1 , m = 1 , . . . , iV-r , n = 1, . . . , N y and p = 1, . . . N z , subject 



0, 



(20) 



SyU 



yu , m,N y ,p 



S z u ri 



0. 



Here, part of the following notation has been employed for the sake of simplicity: 
(21) 
5 t ul 



= u k+1 

"m^ri.p rn,n,p 

X2 u k _ „,k 

X 111,71, p 



8 2 u k 

u z ^'m,n,p 



..k-l X2 k _ k+1 _ k \ y k-l 

m,n,pi t m,n 7 p m,n 7 p m 7 n,p 1 m,n,p> 

ni n 2?/ 4- 77 fi^tt = ?/ 2?/^ -4- 77 

m+l,n,p m,n,p 1 m — l,n,p' y m,n,p m,n+l,p m,n,p 1 m,n — 1 ,p J 



m,n,p+l 



^m,n,p 



U 



m,n,p- 
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tk- 



tk 



tk+1 




Figure 2. Forward-difference stencil depicting implicit finite- 
difference scheme (fT9|) at site (m,n,p) at time tk- The circles 
represent known data at the fcth iteration of the method, while the 
crosses denote the unknown variables. 



The forward-difference stencil of the method is presented in Fig. [2] for convenience. 
Moreover, we introduce the composite operators 5t x = 8 x $t, $ty — 5 y 5t, and Stz — 
S z 8t, and the constant Av — AxAyAz. In these terms, the Hamiltonian of the 
lattice site at position (to, n,p) and the total energy of the system are, respectively, 
(22) 



rrk 

m,n,p 



E k 



1 / 7/^+1 — 

- 1 - I m,n,p rn,n,p 



At 



k \ 

m,n,p/ 



+ 



11V 



(Ax) 2 
r k+i \2 , t k \5 



(Ay) 2 



(Azf 

V(u k+1 ) + V(u k ) u k+1 +u k 

1 J 1 ; 



N 



N 



m,n.p— 1 i t j= 1 

, (<5,<to)(<W ? .o)" ' 



(5 x u^)(5 x u k ) t (SyU^MSyU*) 



(Ax) 



(Ay) 5 



(Azf 



Av. 



In the following and for the sake of simplification, we will set TV = N x 



N„ 



N, 



Proposition 5. The discrete rate of change of energy with respect to time of system 
(JT9]) subject to the discrete boundary conditions (|20|) at the kth instant of time is 
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given by 
(23) 
E k _ E k-i 

At 



5xU Q,i,J 5tU o,id , 5 v u to,j s tUi, ,j t 5 z u k i ] S t ul jfi 



(Ax) 2 2At (Ay) 2 2At (Az) 



2At 




2AxAt 



2AyAt 



dtzUm t n t p-i 

2AzAt 



{^xu^ id ){5tu^ itj ) (Styul 0d )(S t ul 0jj ) (Stzul jfi )(6 t ul ji0 



N 

-i E 

m,n,p— 1 



(2AxAt) 2 

( $t u m ,n,p 

2At 



(2AyAt) 2 



+ 



(2AzAt) 2 



Av. 



Proof. For the sake of simplification, we adopt the convention 5fu k 



t a m,n,p ^m^n^p 



fe+1 



u m,n, P for every k = 0, 1, ... M - 1, m = 0, 1, . . . , N x + 1, n = 0, 1, . . . , N y + 1 
p = 6, 1, . . . , N z + 1. Notice then that 



1 / 5t U m,n,p 

2 \ At 



1 / u t a m. n 



At 



U m,n,p &t U m,n,p 



and 



^( u rrt,n,p) ^( u m,n,p) ^( u m,n,p) ~f" ^( u m,n,p) _ ^( u rn^n,p) G(u mnp ) $t u n 



m,7i,p 



,k+l 



,k-l 



l m,n,p u>m,n,p 



Moreover, 
(*) 

(3xU mn p)(5 x u mnp ) (dxU mn p)(5 x u mnp ) 



2(Ax) 2 



2(Ax) 2 



3x U m,n,p 

Ax 



StxU 



m,n,p 



2Ax 



^x u m,n,p &tU m n p $xU m n p 0~tU m+l n p S x U m _ l n p StU m n ] - 



(Ax) 2 2 (Ax) 2 2 (Ax) 2 2 ' 

which is an expression that may be further simplified for computational purposes 
as a consequence of the convention 



JxU n 



3x u m,n,p ^t u m+l,n,p 



(Ax) 2 2 

Summing the identities Q over all indexes m, n and p, noticing the presence of 
a telescoping sum, and applying the discrete Neumann boundary condition, we 
obtain the following sequence of equalities: 



(^x u nt,n,p)(^x u m ,n,p) {^x'U mnp )(S x U mn ^ p ) 

m,n : p—l 



2(Ax) 2 



2(Ax) 



N 



N 



E^x u m.n,p 3t u m,n,p 
(Ax) 2 2 ^ 

,n,p— 1 n,p—l 

N 



N 



(jx u m,n,p jxU." 



m — l,n,pj 



E 



- (Ax) 2 



fixUQ n p 5tu\ n p 
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In similar fashion, one can verify that substituting the difference 5 X for S y or S z in 
both sides of this last equation yields a valid equality. It is now straight-forward to 
verify that 



N 

E 

m,n,p=l 



At 



N 5<u k 

E u t a m,n,p 
, 2At 

m,n,p—l 



w t Uj m,n,'p 



+ 



(At) 2 

^( u rn,n,p) ~ ^*( M m,n,p) 



,fc+l 



fe-1 



(Ax) 2 

N 

E 



+ 



(Ay) 2 



(Az) 2 



(Ax) 2 2 At 



^Xo.j ^"i,l,3 j,0 <*«"i,j,l 

' (Ay) 2 2At (Az) 2 2Ai 



TV 



1 ^ (2At) 2 

m,n,p— 1 



(Ax) 2 



-7 



N 

E 



<5 t u 



AT 



m,n.p 



2At 



E 



(Ay) 2 (Az) 2 
(Ax) 2 2At 



(Ay) 2 2A< (Az) 2 2At 

Commutativity of the discrete operators yields S t u^ np S t Slu^ np = 8 t u k mnp {5 t u k m+lnp - 2<5 t u ? fe n n p + «m-i,„, p ) 
Summing over all indexes m = 1, 2, . . . , TV, applying the discrete version of Green's 
first identity and summing next over all indexes n,p = 1, 2, . . . , N, it follows that 

JV AT AT 

E 5 * U m,n,p*t<^«m,n,p = ~ E ^ U 0,nj,*t«0,nj, ~ E (^"m-l.nj.) ! 
m,n,p— 1 n,p—l m,n,p— 1 

moreover, similar relations may be obtained for the cases when the operator # x is 
replaced by 5 y or S z . In this circumstances, we readily obtain that 



At 



N 

E 

m,n,p— 1 



2AxAt 



+ 



J ty a m ,n-l,p 

2AyAt 



+ 



2AzAt 



Av 



N 



E 



(**x«o,i,j)(^«o > i,j) . (*t»«i,oj)(*t«?,o,j) . (^ u i,j,o)(^ u i,j,o) 



-7 E 



(2AxAt) 2 

°t u m,n,p 



+ 



{2AyAty 



(2AzAtf 



Ai 



N 



+ 



2At 



Av- £ 



*i3=l 



(Ax) 2 2At (Ay) 2 2At 



(Az) 2 2At 

(^»<,-,o)(^<j,o) 



Aw + ^ 



2(Ax) 2 At 



+ 



Aw, 



2(Az) 2 At 

whence the result follows after an easy simplification in the last two sums. 



2(Ay) 2 At 



□ 
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A direct comparison between this result and Propositions Q] and |4] shows that 
the numerical method presented here consistently approximates the derivative of 
the total energy of a problem described by either (JlJ or ([13"]) . As a consequence, 
the method proposed in this work is capable of preserving the total energy of a 
conservative system described by the continuous equation ([1} or the discrete system 
(|13|) , which is a physical scenario that appears when both /3 and 7 are equal to zero 
and when a fixed boundary is considered. 

For the next result, stability means stability order n (see [37] ). 



Proposition 6. Let V' be identically equal to zero, and let J 
scheme (|19j) to be stable it is necessary that the condition 

1 1 1 



0. In order for 



(24) 



[(At) 2 - 0At] - [7 + m 2 At] At < 4 



{Ax) 2 (Ay) 2 (Az) 2 
be satisfied. 

Proof. For every m,n,p = 1, . . . , N and every k = 0, . . . , M— 1 , let u \ +1 



= u k+1 

m,n,p m,n.p 

and u^t} „ n = ut, „ and let ujL „ „ be the two-dimensional column vector whose 



components are u\ mnp and m „ p 
as 



(u fc+1 ) 



TV 

m,n,p— 1 



In these terms, scheme (fTU|) can be presented 

(Ai)a r « , ^ , « 



(Ax) 2 



(Ay) 2 



(Az) 2 



U 







™,n,PJ m,n,p=l ' 



where 



1 - 



(3At 
f3At 



si 




si 


-yAt 


m 2 (At) 2 


(Ax) 2 


(Ay) 2 


(Az) 2 _ 


+ — " 


- 2 




S l 


sl 


7At 


m 2 (At) 2 


(Ax) 2 


(Ay) 2 


(Az) 2 




- 2 



We apply Fourier transform in order to reach the expression 

7 ^(i- 2 ^) 2 (-r 2i 



fe+i _ 



1 



(Az) 




where the 'hat' operator obviously denotes Fourier transform. We identify the 
2x2 matrix multiplying n in the above equation as the amplification matrix 
^4(£> C> •>) °f our problem. Moreover, it is easy to check that the eigenvalues of this 
matrix when £, £ and ? are all equal to 7r, are given by 



1 - 2R 2 (At) 2 ± J (I - 2R 2 (At) 2 ) 2 - /i(7r,7r,7r)5(7r,7r,7r) 



5(7r,7r,7r) 



where 



1 



(Ax) 2 (Ayf 



Suppose for a moment that 1 — 2R 2 (At) 2 < 



(Az) 2 ' 

-g(-rr,ir,ir) 



If the radical in the 

expression above yields a pure real number then |A_| > 1. So for every positive 
integer I, \\A l \\ > grows faster than K\ + IK2 for any constants K\ and Ka- 

A similar situation prevails when the radical is a pure imaginary number, except 
that in this case | • | represents the usual Euclidean norm in the field of complex 
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Figure 3. Graph of the approximate solution U60,60,60 versus 
time, obtained by driving system (fl3|) at a frequency equal to 0.9 in 
the forbidden band-gap of the continuous-limit medium, and two 
different amplitudes: 1.42 (left) and 1.43 (right). All other param- 
eters are set equal to zero, a time step equal to 0.05 was fixed, and 
a system with an absorbing boundary consisting of 200 3 coupled 
nodes was considered. The graphs are presented as evidence of the 
presence of supratransmission in the system under study. 



numbers. Therefore in order for our numerical method to be stable it is necessary 
that 1 — 2R 2 (At) 2 > —g(n, tt, tt), which is what we wished to establish. □ 



Corollary 7. Let V' be identically equal to zero, let J = 0, and suppose that 
Ax — Ay — Az. In order for scheme (|19j) to be stable it is necessary that the 
condition 

(25) (12i? 2 - m 2 ){At) 2 - (7 + \2(3R 2 )At < 4 

be satisfied, for R = 1/ Ax. □ 

It is important to notice that the order of consistency of the method as defined 
by the truncation error is 0((Ax) 2 + (Ay) 2 + (Az) 2 + (At) 2 ). Also, it is worth 
mentioning that we have approximated V (u(mAx, nAy,pAz, kAt)) in the right- 
hand side of Eq. (TTJ) through the discrete derivative of V with respect to u presented 
in (|19p and not through the direct evaluation V'(u^ nnp ) in view of the fact that 
such standard scheme is known to be highly unstable [55] , 



3.2. Spherical problem. In order to approximate solutions to © in the closure 
of the solid sphere with center in the origin and radius L, we let e = T"o < ri < • • • < 
tm+i = L and — to < ti < • • ■ < tjsr — T be regular partitions of [e, L] and [0, T], 
respectively, into M + 1 and N subintervals of lengths Ar and At, respectively. 
Denote the approximate value of v(rj,tk) by Vj. The finite-difference scheme used 
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to pursue that task is the implicit method 

(At) 2 (Ar) 2 7 2At 

fc-n 



2At (Ar) 2 2 J 



(e + jAr) m ^ J(e + jAr) = 0, 

defined for every j = 1, . . . , M and every k = 1, . . . , N — 1. Meanwhile, the total 
energy of the system (1101) and the Neumann boundary condition of the problem 
(|T2"j) at the fc-th time step will be respectively approximated using the schemes 
(27) 




and 

/ 9 o\ ^M+l ~ , "M+l + V M _ n 

1 j Ar + 2(e + MAr) 2 ' 

The following results summarize the most important numerical properties of the 
method proposed in the present subsection. The proofs are omitted in view that 
they are similar to the corresponding results of the Cartesian case. 

Proposition 8 (Macfas-Diaz et al |39j ) . The discrete rate of change of energy with 
respect to time of system (|26|) subject to the discrete boundary condition (|28[) at the 
kth instant of time is given by 
(29) 

At 2f^l 2At [ At{A r y 




Ar 



for e equal to zero. □ 

Proposition 9 (Maci'as-Dfaz et al [39]). Let V be identically equal to zero, and let 
J = 0. In order for finite- difference scheme @ to be stable order n it is necessary 
that 

f At\ 2 At At ? (At) 2 

30 — <l + 7 — +13 _+ m 3i_2_ 

\ArJ 4 (Ar) 2 4 

□ 
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Driving Amplitude 

FIGURE 4. Bifurcation diagram of total energy of node located at 
sine (60, 60, 60) over a time period of 200 versus driving amplitude, 
for a diving frequency equal to 0.9 in the forbidden band-gap region 
of the continuous-limit medium. All other parameters are equal 
to zero, a time step equal to 0.05 was employed, and a system 
with an absorbing boundary consisting of 200 3 coupled nodes was 
considered. The graph is presented as evidence of the presence of 
supratransmission in the system under study. 

We must remark that the computational technique presented in this section 
makes use of Newton's method to approximate solutions of systems of nonlinear 
equations. In each iteration, the system of equations derived from Newton's method 
is linear, and may be solved using Crout's technique for tridiagonal systems. 

4. Applications 

Throughout this section, we suppose that the function driving the boundary as- 
sumes the expression <fi(t) = Asin(ftt), where the driving frequency takes on values 
in the forbidden band-gap region of the continuous-limit medium < y/m? + 1. 
The explicit scheme resulting from setting /3 equal to zero in (|19|) is employed to 
verify the validity of our results in the case of weak internal damping. 

4.1. Nonlinear supratransmission. The process of supratransmission in nonlin- 
ear systems submitted to harmonic driving is completely characterized by a sudden 
increase in the energy injected into the system by the driving source [Tj. Thus, 
the method employed to determine the critical value at which supratransmission 
starts, given a fixed frequency ft in the forbidden band-gap of the system, consists 
in computing the associated energy E for various driving amplitudes A in an inter- 
val containing the supratransmission threshold, over a fixed, relatively long period 
of time; in these circumstances, the graph of E versus A will evidence a point of 
discontinuity where a drastic increase in the total energy of the system takes place. 

In the case of an infinite number of coupled junctions (u m .n,p)m n p=i satisfying 
problem (p~3|) . we consider a large finite subsystem consisting of N x N x N cou- 
pled junctions, with damping coefficient 7m,n,p including the effect of an absorbing 
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boundary in the farthest junctions from the point of intersection of the three driving 
boundaries. That is, we let 



where 1 < N Q < N. 

Let all constant parameters in the differential equations of (ITU)) be set equal 
to zero. Following the method described in the previous paragraph, we submit 
system (|13[) to harmonic driving with frequency in the forbidden band-gap of the 
continuous-limit medium over a time period of [0, 200], and compute the associated 
total energy of the system for several driving amplitudes. In order to avoid the 
generation of shock waves at the origin around the time t — 0, we opt for slowly and 
linearly increase the driving amplitude from to its actual value A. Numerically, 
we choose a time step equal to 0.05 (so that the stability condition provided by 
Proposition [B] is clearly satisfied), fix a cubic system of dimension 200 3 , and we 
damp the farthest nodes from the origin using (|3"Tj) and No = 50. 

To start with, fix a driving frequency of 0.9 and chose two different amplitude 
values: A = 1.42 and A = 1.43. The time behavior of the solution of the node 
located at site (60, 60, 60) as a result of driving system (Q~3|) under the circumstances 
described in the paragraph above for the two amplitudes considered is displayed in 
Fig. [3] It is worth noticing that the wave signals transmitted into the system for a 
driving amplitude equal to 1.42 posses a very low amplitude when compared against 
the driving amplitude itself (left graph). On the other hand, a driving amplitude of 
1.43 produces wave signals of higher amplitude (right graph). Moreover, the total 
energy at site (60, 60, 60) is computed by integrating the discrete Hamiltonian, 
obtaining, in the first case, an energy equal to 68.7613, while a total energy of 
161.3648 is obtained in the second, whence the existence of a critical amplitude 
between 1.42 and 1.43 at which supratransmission starts is suspected. 

Next, we compute the total energy at site (60, 60, 60) of system (|13| for several 
amplitude values around the suspected critical values, and for a fixed frequency 
equal to 0.9. As before, all other parameters are set equal to zero, and we use the 
same numerical setting as before. In these circumstances, we present the graph of 
total energy versus amplitude on the time period [0, 200]. The results are presented 
in Fig. |4] and confirm that a drastic increase in the total energy of the node 
appears for a driving amplitude between 1.42 and 1.43, proving thus the presence 
of nonlinear supratrasmission in our system, at least for a driving frequency of 0.9. 
In this point it must be mentioned that we have established that supratransmission 
is likewise present for other choices we made of the parameter Q. 

4.2. A counter-example. Numerically, let us fix a time step At = 0.02, and let 
the radial step Ar and the parameter e both equal At. In all our computations we 
consider a fixed time period. Moreover, in order to avoid the generation of shock 
waves, the driving amplitude will increase linearly and slowly from to its actual 
value A during a relatively short period of time, before the initial instant t = 
takes place. 

Assume that the medium has no damping. In order to simulate an unbounded 
medium, we approximate solutions to problem ((9|) in a closed sphere S with center 



(31) 
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Figure 5. Total energy versus driving frequency fi and driving 
amplitude A in undamped system @ with J = and potential 
equal to that of a nonlinear Klein-Gordon equation (left), and equal 
to that of a classical sine- Gordon medium (right). The harmonic 
driving function at the origin is defined by <fi(t) — Asm(£li), and a 
time period of 200 was fixed. 

in the origin and radius L — 6, in which the parameter 7 slowly increases in 
magnitude from to 1 outside the open sphere with center in the origin and radius 
5, simulating thus an absorbing boundary. More precisely, we let 



It is worth noticing that if the potential function V is that for a sine-Gordon sys- 
tem, then rV'(v/r) is approximately equal to zero for nonzero values of r sufficiently 
close to zero. It is therefore expected that the medium behaves in a linear fashion 
around the origin and, particularly, that the medium does not support the process 
of supratransmission under the presence of harmonic perturbations at the origin. 
Of course, this claim will be confirmed numerically next for both Klein-Gordon and 
sine-Gordon systems. 

Throughout, we let A range in [0, 20]. Let us take a fixed frequency f2 = 0.9 in 
the forbidden band-gap of a continuous Klein-Gordon medium described by model 
(jH) . The associated total energy of the system during the fixed period of time is 
computed, obtaining evidence of a continuous increase in the amount of energy 
injected in the system by the driving boundary, with no apparent discontinuities 
in the total energy of the system. Next, we let f2 take on values in the interval 
[0,1]. The graph of total energy versus driving amplitude and driving frequency 
is presented as Fig. [5Ja), and the results evidence that the total energy increases 
smoothly as the driving amplitude is increased. This fact supports our claim that 
the process of nonlinear supratransmission is not present in this medium. 

We proceed to examine the sine-Gordon case next. Preliminary computational 
results show that, for a fixed frequency il = 0.9 in the forbidden band-gap, the total 
energy of the system during a fixed period of time equal to 20 increases smoothly 
with respect to A. Thus, we let f2 range between and 1 and, for each pair (fi, A), 
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Figure 6. Time evolution of a localized solution of (fTB")) with 
7 = 0.005, (3 = m 2 = and J = 0.01, subject to harmonic 
driving with a frequency f2 = 0.9 in the forbidden band-gap of 
the continuous-limit system. The driving function takes the form 
(f>(t) = A(t) sin(f2i) where A is given by (f34|) . and bi = 1 for every 
i = 1, 2, 3, 4. The graphs show the time evolution of regions of high 
energy. 
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Time 



Figure 7. Time-dependent graph of local energy density H of 
the node located at site (60, 60, 60) in a medium described by fp~3|) 
with 7 = 0.005, j3 = m 2 = and J = 0.01, as a consequence 
of driving it harmonically on the boundary at a frequency of 0.9 
in the forbidden band-gap of the continuous-limit medium, and a 
driving amplitude (l34l) with b\ = b% = b$ = &4 = 1. The peaks 
of the graph record the propagation of localized nonlinear modes 
produced by the four nonzero bits transmitted into the medium. 

compute the associated total energy of the system. In this context, Fig. E](b) 
prescribes the total energy of the system versus f2 and A. Our results show that, 
the phenomenon of nonlinear supratransmission is absent in the case of radially 
symmetric sine- Gordon systems. 

4.3. Propagation of signals. The study of localized nonlinear modes in (1 + 1)- 
dimensional sine-Gordon systems is a topic of research that has produced a large 
amount of valuable results. Nowadays, the specialized literature in the field pos- 
sesses results on this topic that range from the analytical aspects of the problem, 
to the numerical, to the physical, including those works where the propagation of 
localized modes are studied in relation with the process of nonlinear supratrans- 
mission El El HQ]- 

In our study, it is particularly important to recall that supratransmission in 
semi-unbounded, sine-Gordon systems subject to Dirichlet harmonic driving has 
been characterized by the generation of moving breathers at the driving boundary 
once the driving amplitude has reached its critical value, and a method to control 
the propagation of these modes in such systems has been proposed [30] . Thus, it 
seems natural to generalize this technique to the case of (3 + l)-dimensional, semi- 
unbounded systems governed by sine-Gordon equations and subject to harmonic 
driving at the boundary. 

Following the method proposed in [40], let us fix < 1, and let <j>(t) = A(t) sm{Q£). 
We let P be a multiple of the driving period, assume that a bit b € {0, 1} will be 
transmitted into a medium (|13[) during the period of time [0,P], and suppose that 
A a represents the critical amplitude at which supratransmission starts for the fre- 
quency f2. Define the driving amplitude function Aq as 

(33) A (t; b) = Cb (e~ m / 2 5 - 6 -™/oa 5 ^ > 
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where C is a positive real number depending on VL that works as an amplification 
factor. The idea behind the definition of A® is that for the process of nonlinear 
supratransmission requires a certain amount of time to start to irradiate energy 
into the medium, during which the driving amplitude must take on values above 
the critical value A s . 

Let XA(t) represent the characteristic function on the set AC1 evaluated at t, 
which is equal to 1 if t G A, and is equal to zero otherwise. In our study, we will 
fix 7 = 0.005, (3 = m 2 = and J = 0.01, for which the value of amplitude at which 
supratransmission starts is A s = 1.41. Numerically, we choose a step size for time 
equal to 0.005, and fix a bounded cube of sides equal to 200. Moreover, the period 
P of signal generation will be equal to 150, the amplification factor will be equal 
to 3, and the binary sequence (b\, . . . , bk) will be transmitted into the medium by 
means of the harmonic driving function <p(t) with amplitude function defined by 

k 

(34) A(t) =J2Mt-(i-l)P;bi)X[(i-i)P,iP]{t)- 

i=l 

The binary sequence (1,1,1,1) will be transmitted into (|13|) using the amplitude 
function just defined. In these conditions, Fig. [5] presents the time evolution of 
the local energy density of a localized solution of the medium as a result of being 
subject to harmonic driving with amplitude defined by (|34[) . during the first period 
of generation of signals. The blue zones presented in these graphs are regions of 
high energy produced at the origin. As time evolves, these regions clearly expand 
and move away from the origin around the line x — y ~ z = t, for t G R + . 

Finally, Fig. [7] presents the time behavior of the local energy density at site 
(60,60,60) in system (fl3|) . as a result of transmitting the binary code (1,1,1,1) 
by means of perturbations on the boundaries. It is worth noticing that each of 
the peaks in the graph is a result of the localized traveling solutions — produced 
by each of the nonzero bits generated at the boundaries — , which passes by the 
node at site (60,60,60) and moves away from the origin. Moreover, there exists a 
gap in time approximately equal to 150 between two consecutive peaks, which is in 
perfect agreement with the value of the period P of signal generation. The results 
evidence the possibility of accurately transmitting binary information into system 
(JXJ) through suitable perturbations of the driving boundary. 

It must be mentioned that similar results (not included here) are obtained for 
the approximation to the continuous case described by (p}, proving thus that the 
presence of nonlinear supratransmission in (3 + l)-dimensional, dissipative sine- 
Gordon equations does not depend on discreteness. 

5. Conclusion 

In this work, we have presented conditionally stable, finite-difference schemes 
that consistently approximate the solutions to problems (p}, ((9J and fp~3|) . Asso- 
ciated with these schemes, we have introduced discrete schemes to approximate 
consistently the local energy densities of the media and their total energy func- 
tions, in such way that the corresponding discrete rates of change of energy with 
respect to time consistently approximate their respective continuous rates of change 
of energy. In particular, if no dissipation is present and under suitable boundary 
conditions, the proposed methods are conservative. 
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Also, we have provided relevant numerical evidence that the process of nonlinear 
supratransmission is not present in media described by undamped radially sym- 
metric sine-Gordon equations perturbed harmonically at the origin, proving thus 
that not every nonlinear system with a forbidden band-gap for the frequency in 
the linear dispersion relation is able to sustain this nonlinear process (contrary to a 
conjecture in the literature pQ). Our computations are supported empirically in the 
case of a sine- Gordon system by the fact the differential equation in ([9]) is approx- 
imately linear close to the origin, and analytically by the well-known fact that the 
origin of such systems is incapable of creating localized coherent structures. On the 
other hand, a similar three-dimensional medium (discrete or continuous), bounded 
in the first octant by the coordinate planes and subject to harmonic driving of the 
Dirichlet type on the boundaries, does exhibit supratransmission. Our results (pre- 
sented in Section 14. 1| show a well-defined occurrence of the critical value at which 
supratransmission starts. 

It is interesting to notice that supratransmission in discrete 1-dimensional chains 
of oscillators in achieved when the first oscillator is harmonically perturbed at a fre- 
quency in the forbidden band-gap. In a 2-dimensional scenario, supratransmission 
is achieved when the boundary lines of a semi-unbounded domain are perturbed at 
the right frequency [BJ. Similarly, a semi- unbounded region in the 3-dimensional 
case presents supratransmission when the boundary surfaces are perturbed at fre- 
quencies in the forbidden band-gap. Following this pattern, an (n+ l)-dimensional 
semi-unbounded system of oscillators described by coupled sine-Gordon equations 
may present supratansmission when the n+1, ?i-dimensional boundaries are subject 
to harmonic driving with a frequency in the forbidden band-gap. 

Another application to the transmission of localized nonlinear modes in (3 + 1)- 
dimensional systems governed by continuous sine- Gordon equations was provided 
in this work. The system was the same semi-unbounded medium studied before — 
the medium governed by (fTJ), defined in the first octant and driven harmonically at 
the coordinate planes. By making use of nonlinear supratransmission, our results 
(summarized in Section [4.3l) show that a controlled propagation of wave signals can 
be achieved. Moreover, the propagating nodes are seen to be traveling breathers 
that move away from the origin on the line x = y = z = t, for t G R + , which is in 
perfect agreement with the (1 + l)-dimensional scenario [40] . 
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Appendix A. Computational setting 

There are several important remarks on finite-difference scheme (|19|) associated 
with partial differential equation (|TJ) . For the sake of simplification, we will assume 
that the spatial step-sizes Ax, Ay and Az are equal, and that N x , N y and N z are 
all equal to N. Moreover, we will let J be equal to zero. 

• First of all, if /? is equal to zero and V is identically equal to zero then 
the resulting differential equation is the damped, linear Klein-Gordon-like 
equation. In such case, the finite-difference method obtained is explicit. 

• If f) is equal to zero but V is not the function identically equal to zero, 
the partial differential equation obtained is a damped, nonlinear Klein- 
Gordon-like equation. In this case, method (|I9j) is nonlinear and explicit; 
in fact, an application of Newton's method is needed to obtain the value 
u rnn p & om the known approximations at times k — 1 and k, for every 
m,n,p — 1, 2, . . . , N (this was the case when performing the application in 
Section I3~I]>. 

• Let (3 be a positive real number. If V' is identically equal to zero then the 
resulting partial differential equation is a damped, linear equation, while 
the finite-difference scheme associated with it is likewise linear and implicit. 
Meanwhile, if V is not identically equal to zero then an application of 
Newton's method is indispensable; we will describe this last scenario in 
more detail now. 

Notice first of all that the algorithm of division implies that, for ev- 
ery positive integer i = 1,2,..., (N + 2) 3 , there exist unique nonnegative 
integers m, n and p such that 

(35) i = m(N + 2) 2 + n(N + 2) + p + 1. 

Conversely, if to, n,p € {0,1,...,JV + 1} then the value of i given by the 
formula above is in the set {I, 2, . . . , (N + 2) 3 }. Thus the term may 
be unambiguously represented by y\, and the left-hand side of Eq. (fT9|) 
will be denoted by //\ 

Following Newton's method, the approximations at times k — I and k 
are assumed to be known in order to compute the approximation at the 
(k + l)st time. The Jacobian of the problem is a sparse matrix; indeed, 
notice that for each triplet (m, n, p) with m,n,p = I, 2, . . . , N, the following 
are the only nonzero partial derivatives: 

(36) 

df[ _ dft _ dtf _ df[ _ df[ Off (3 

dyi+i %-i dy l+N+2 dyi_ {N+2 ) dy l+{N+2) 2 dy t _ {N+2 y 2At(Axf 

and 

(V) d A( \ = J_ 3 1 7 . m 2 {z-u k i - 1 )V'{z) + V{u k r 1 )-V{z) 

[ ' d yi [ > {Aty + At{AxY + 2At + 2 + (z-ut 1 ) 2 
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